0%

【工程优化】最优化算法--牛顿法、阻尼牛顿法及单纯形法

发表于 更新于
本文字数： 5.8k 阅读时长 ≈ 8 分钟

牛顿法

使用条件：目标函数具有二阶导数，且海塞矩阵正定。

优缺点： 收敛速度快、计算量大、很依赖初始点的选择。

算法的基本步骤：

已知目标函数$f(x)$，梯度$g(x)$, Hessan矩阵$G(x)$, 给定误差限$\epsilon$:

  • 步骤1：选定初始点$x_0$，计算$f_0=f(x_0), k=0$;
  • 步骤2：计算$g_k=g(x_k)$，如果$\Vert g_k\Vert\le\epsilon$, 算法停止, $x^\star=x_k$，否则转到步骤3：
  • 步骤3：计算$G_k=G(x_k)$, 由方程$G_kd^k=-g_k$, 解得$d^k$;
  • 步骤4：令$x_{k+1}=x_k+d^{k}, k=k+1$, 转到步骤2

算法流程图：

​ 图1

阻尼牛顿法

与牛顿法基本相同，只是加入了一维精确搜索。在牛顿迭代中，取$d^{(k)}=-[\nabla^2f(x^{(k)}) ]^{-1}\nabla f(x^{(k)}) $ , 加入精确一维搜索：$\min f(x^{(k)}+\lambda_kd^{(k)})$, 求得$\lambda_k$, 然后更新：$x^{(k+1)}=x^{(k)}+\lambda_kd^{(k)}$.

优缺点：改善了局部收敛性。

我们假设要求$f=(x-1)\cdot(x-1)+y\cdot y$的最小值，具体算法实现如下，只需要运行NTTest.m文件，其它函数文件放在同一目录下即可：

1、脚本文件NTTest.m

1
2
3
4
5
6
7
8
9
10
11
12
13
clear all
clc
 
syms x y
f=(x-1)*(x-1)+y*y;
var=[x y];
x0=[1 1];eps=0.000001;
 
disp('牛顿法:')
minNT(f,x0,var,eps)
 
disp('阻尼牛顿法:')
minMNT(f,x0,var,eps)

2、minNT.m

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
 function [x,minf]=minNT(f,x0,var,eps)
%目标函数：f
%初始点：x0
%自变量向量：var
%精度：eps
%目标函数取最小值时的自变量值：x;
%目标函数的最小值：minf
 
format long;
if nargin==3
    eps=1.0e-6;
end
tol=1;
syms L
% x0=transpose(x0);
while tol>eps %不满足精度要求          
    gradf=jacobian(f,var);      %梯度方向
    jacf=jacobian(gradf,var);   %雅克比矩阵
    v=Funval(gradf,var,x0);%梯度的数值解
    tol=norm(v);%计算梯度（即一阶导）的大小
    pv=Funval(jacf,var,x0);%二阶导的数值解
    p=-inv(pv)*transpose(v);    %搜索方向
    x1=x0+p';%进行迭代
    x0=x1;
end
 
x=x1;
minf=Funval(f,var,x);
format short;

3、minMNT.m

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
function [x,minf]=minMNT(f,x0,var,eps)
%目标函数：f
%初始点：x0
%自变量向量：var
%精度：eps
%目标函数取最小值时的自变量值：x;
%目标函数的最小值：minf
 
format long;
if nargin==3
    eps=1.0e-6;
end
tol=1;
syms L
% x0=transpose(x0);
while tol>eps %不满足精度要求          
    gradf=jacobian(f,var);      %梯度方向
    jacf=jacobian(gradf,var);   %雅克比矩阵
    v=Funval(gradf,var,x0);%梯度的数值解
    tol=norm(v);%计算梯度（即一阶导）的大小
    pv=Funval(jacf,var,x0);%二阶导的数值解
    p=-inv(pv)*transpose(v);    %搜索方向
    %%%%寻找最佳步长%%%
    y=x0+L*p';
    yf=Funval(f,var,y);
    [a,b]=minJT(yf,0,0.1);
    xm=minHJ(yf,a,b);           %黄金分割法进行一维搜索最佳步长
    x1=x0+xm*p';%进行迭代
    x0=x1;
 
end
 
x=double(x1);
minf=double(Funval(f,var,x));
format short;

4、minHJ.m

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
function [x,minf]=minHJ(f,a,b,eps)
%目标函数：f
%极值区间左端点：a
%极值区间右端点：b
%精度：eps
%目标函数取最小值时自变量的值：x
%目标函数所取的最小值：minf
 
format long;
if nargin==3
    eps=1.0e-6;
end
 
l=a+0.382*(b-a);            %试探点
u=a+0.618*(b-a);            %试探点
k=1;
tol=b-a;
 
while tol>eps&&k<100000
    fl=subs(f,findsym(f),l);        %试探点函数值
    fu=subs(f,findsym(f),u);        %试探点函数值
    if fl>fu
        a=1;                        %改变区间左端点
        l=u;
        u=a+0.618*(b-a);            %缩短搜索区间
    else
        b=u;                        %改变区间右端点
        u=l;
        l=a+0.382*(b-a);             %缩短搜索区间
    end
    k=k+1;
    tol=abs(b-a);
end
if k==100000
    disp('找不到最小值！');
    x=NaN;
    minf=NaN;
    return;
end
x=(a+b)/2;
minf=subs(f,findsym(f),x);
format short;

5、minJT.m

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
function [minx,maxx]=minJT(f,x0,h0,eps)
%目标函数：f
%初始点：x0
%初始步长：h0
%精度：eps
%目标函数取包含极值的区间左端点：minx
%目标函数取包含极值的区间右端点：maxx
 
format long
if nargin==3
    eps=1.0e-6;
end
 
x1=x0;
k=0;
h=h0;
while 1
    x4=x1+h;        %试探步
    k=k+1;
    f4=subs(f,findsym(f),x4);
    f1=subs(f,findsym(f),x1);
    if f4<f1
        x2=x1;
        x1=x4;
        f2=f1;
        f1=f4;
        h=2*h;      %加大步长
    else
        if k==1
            h=-h;   %方向搜索
            x2=x4;
            f2=f4;
        else
            x3=x2;
            x2=x1;
            x1=x4;
            break;
        end
    end
end
 
minx=min(x1,x3);
maxx=x1+x3-minx;
format short;

6、Funval.m

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
function fv=Funval(f,varvec,varval)
var=findsym(f);
varc=findsym(varvec);
s1=length(var);
s2=length(varc);
m=floor((s1-1)/3+1);
varv=zeros(1,m);
if s1~=s2
    for i=0:((s1-1)/3)
        k=findstr(varc,var(3*i+1));
        index=(k-1)/3;
        varv(i+1)=varval(index+1);
    end
    fv=subs(f,var,varv);
else
    fv=subs(f,varvec,varval);
end

运行结果如下图：

图2

单纯形法

单纯形法的理论还有点复杂，而本文主要针对算法的基本实现，因此，理论部分就此略过，详情可以参考网上的相关资料。下面给出具体的实现：

我们以具体实例来说明：

假定线性规划问题如下：

化为标准型可得：

具体的Matlab实现如下：

1、脚本文件：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
clear all
clc
% A=[2 2 1 0 0 0
%    1 2 0 1 0 0
%    4 0 0 0 1 0
%    0 4 0 0 0 1];
% c=[-2 -3 0 0 0 0];
% b=[12 8 16 12]';
% baseVector=[3 4 5 6];
 
A=[1 1 -2 1 0 0
   2 -1 4 0 1 0
   -1 2 -4 0 0 1];
c=[1 -2 1 0 0 0];
b=[12 8 4]';
baseVector=[4 5 6];
 
[x y]=ModifSimpleMthd(A,c,b,baseVector)

2、ModifSimpleMthd.m文件

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
function [x,minf]=ModifSimpleMthd(A,c,b,baseVector)
%约束矩阵：A
%目标函数系数向量：c
%约束右端向量：b
%初始基向量：baseVector
%目标函数取最小值时的自变量值：x
%目标函数的最小值：minf
 
 
sz=size(A);
nVia=sz(2);
n=sz(1);
xx=1:nVia;
nobase=zeros(1,1);
m=1;
 
if c>=0
    vr=find(c~=0,1,'last');
    rgv=inv(A(:,(nVia-n+1):nVia))*b;
    if rgv>=0
        x=zeros(1,vr);
        minf=0;
    else
        disp('不存在最优解');
        x=NaN;
        minf=NaN;
        return;
    end
end
 
for i=1:nVia            %获取非基变量下标
    if(isempty(find(baseVector==xx(i),1)))
        nobase(m)=i;
        m=m+1;
    else
        ;
    end
end
 
bCon=1;
M=0;
B=A(:,baseVector);
invB=inv(B);
 
while bCon
    nB=A(:,nobase);         %非基变量矩阵
    ncb=c(nobase);          %非基变量系数
    B=A(:,baseVector);      %基变量矩阵
    cb=c(baseVector);       %基变量系数
    xb=invB*b;
    f=cb*xb;
    w=cb*invB;
    
    for i=1:length(nobase)  %判别
        sigma(i)=w*nB(:,i)-ncb(i);
    end
    [maxs,ind]=max(sigma);  %ind为进基变量下标
    if maxs<=0              %最大值小于零，输出解最优
        minf=cb*xb;
        vr=find(c~=0,1,'last');
        for l=1:vr
            ele=find(baseVector==l,1);
            if(isempty(ele))
                x(l)=0;
            else
                x(l)=xb(ele);
            end
        end
        bCon=0;
    else
        y=inv(B)*A(:,nobase(ind));
        if y<=0             %不存在最优解
            disp('不存在最优解！');
            x=NaN;
            minf=NaN;
            return;
        else
            minb=inf;
            chagB=0;
            for j=1:length(y)
                if y(j)>0
                    bz=xb(j)/y(j);
                    if bz<minb
                        minb=bz;
                        chagB=j;
                    end
                end
            end                     %chagB为基变量下标
            tmp=baseVector(chagB);  %更新基矩阵和非基矩阵
            baseVector(chagB)=nobase(ind);
            nobase(ind)=tmp;
            
            for j=1:chagB-1         %基变量矩阵的逆矩阵变换
                if y(j)~=0
                    invB(j,:)=invB(j,:)-invB(chagB,:)*y(j)/y(chagB);
                end
            end
            for j=chagB+1:length(y)
                if y(j)~=0
                    invB(j,:)=invB(j,:)-invB(chagB,:)*y(j)/y(chagB);
                end
            end
            invB(chagB,:)=invB(chagB,:)/y(chagB);
            
        end
    end
    M=M+1;
    if(M==1000000)               %迭代步数限制
        disp('找不到最优解！');
        x=NaN;
        minf=NaN;
        return;
    end
end

运行结果如下图：

图3

关于最优化的更多算法实现，请访问http://download.csdn.net/detail/tengweitw/8434549，里面有每个算法的索引说明，当然也包括上述算法。

坚持原创技术分享，您的支持将鼓励我继续创作！